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Abstract. The basic form of drifting sub-pulses is that of a periodicity whose phase depends (approximately 
linearly) on both pulse longitude and pulse number. As such, we argue that the two-dimensional Fourier transform 
of the longitude-time data (called the Two-Dimensional Fluctuation Spectrum; 2DFS) presents an ideal basis for 
studies of this phenomenon. We examine the 2DFS of a pulsar signal synthesized using the parameters of an 
empirical model for sub-pulse behaviour. We show that the transform concentrates the modulation power to 
a relatively small area of phase space in the region corresponding to the characteristic frequency of sub-pulses 
in longitude and pulse number. This property enables the detection of the presence and parameters of drifting 
sub-pulses with great sensitivity even in data where the noise level far exceeds the instantaneous flux density of 
individual pulses. The amplitude of drifting sub-pulses is modulated in time by scintillation and pulse nulling 
and in longitude by the rotating viewing geometry (with an envelope similar to that of the mean pulse profile). 
In addition, sub-pulse phase as a function of longitude and pulse number can differ from that of a sinusoid 
due to variations in the drift rate (often associated with nulling) and through the varying rate of traverse of 
magnetic azimuth afforded by the sight line. These deviations from uniform sub-pulse drift manifest in the 2DFS 
as broadening of the otherwise delta-function response of a uniform sinusoid. We show how these phase and 
amplitude variations can be extracted from the complex spectrum. 
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1. Introduction 

Very soon after the discovery of pulsars, it was noticed 
that the intensity data sometimes showed "second period- 
icities" within each pulse, and that the relative phase of 
this periodic signal drifted in an organised manner from 



pulse to pulse (Drake & Craft 1968). The striking pat- 



terns made by drifting sub-pulses in two-dimensional pulse 
longitude-time diagrams seemed to be saying something 
about the fundamentals of radio pulsar emission, but the 
question was, and still is, 'what?' Numerous studies of the 
phenomenon were conducted in the first few years of pul- 



sar science (e.g. Taylor ct al. 1969; 


Sutton ct al. 1970|; 


Cole 


197C 


; 
; 


Backer 1970a 




Taylor & Hugucnin 1971; 


Backer 


1973 


Page 1972), their frequency steadily decreasing as 



the limitations of existing online (recording) and offline 
(analysis) equipment were reached. Along with attempts 
at fitting the data with sub-pulses and tracking their drift, 
the Fourier technique developed by Backer in the refer- 
ences above (now known as the longitude-resolved fluctu- 
ation spectrum, or LRFS) came to be perhaps the most 
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widely used means of examining the average properties of 
sub-pulse drift. 

With the availability of modern digital back-ends and 
convenient offline computing power, we believe that there 
is now potential for renewed progress in the field of drifting 
sub-pulses. We describe in this paper a means of analysis 
of drifting sub-pulse data using what may be considered 
an extension of the Fourier work of Backer. The technique 
was first conceived as a means of detecting and char- 
acterising drifting sub-pulses in the large population of 
moderately weak pulsars that have been discovered since 
the 1970s when most studies took place. It is also use- 
ful, when sufficient signal-to-noise ratio is available, for 
the investigation of what might be considered the fine de- 
tails of the phenomenon: deviations from constant am- 
plitude sub-pulses with purely uniform phase evolution. 
We begin with a mathematical description of a drifting 
sub-pulse signal as a function of pulse longitude and pulse 
number (Sect, p) and examine its two-dimensional Fourier 
response (Sect. |). We then describe a technique for mea- 
suring the deviations of a given signal from a uniform pure 
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sinusoid (Sect. 4.1), examine what is expected under the 
usual model of drifting sub-pulses (Sect. 4.2), and use this 
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model as a basis for simulations testing the applicability 
of the methods of analysis described here (Sect. 4.2). 



2. Mathematical Model of Drifting Sub-Pulse 
Signals 

2.1. The Form Of Drifting Sub-Pulses 

The characteristic form of drifting sub-pulses is difficult 
to miss in a suitable plot of high-quality data. Plotted 
as intensity as a function of pulse longitude and pulse 
number, a distinctive "tire tread" pattern is made as the 
sub pulses drift monotonically in longitude with each suc- 
cessive pulse. The form is essentially that of a windowed 
two-dimensional sinusoid. The characteristic period of the 
sub-pulses in pulse longitude (P2) and pulse number (P3 
when expressed as a time intervaj^ may be measured by 
a variety of techniques and is usually the first step in any 
analysis of pulsars exhibiting drifting sub-pulses. 

Naturally the actual signal is not a uniform, pure two- 
dimensional sinusoid. The intensity at a given instant de- 
pends on the pulse longitude (the overall emission is, af- 
ter all pulsed at the spin period Pi) and also has a time 
variability that may be intrinsic to the pulsar or induced 
by interstellar scintillation. The spacing of the sub-pulses 
may also vary in a systematic manner with time and/or 



pulse longitude. We shall show later (Sect. 4.2) that un- 
der the most popular model for the origin of drifting 
sub-pulses, we may to a very good approximation model 
these variations in sub-pulse spacing as as a dependence 
of P2 on pulse longitude and/or a time- variability of P3. 
Observational results showing that (in specific pulsars) 
P2 does not vary 



et al. 



197S; 


Lync & Ashworth 1983 


) and that P3 is indcpcn- 


dent of pulse longitude (Backer 1970a|Jb|; 


Krishnamohan 


198C; Wright 1981; Davics et al. 1984 


; Biggs et al. 1987) 



serve to confirm this assertion. The only reported excep- 
tion, to our knowledge, is the sub-pulse modulation re- 
ported in the putative pulsar remnant of Supernova 1987A 
( Middlcditch et al. 200Cl| ), which, unlike other pulsars with 
drifting sub-pulses, is probably caused by precession. We 
therefore proceed to describe the signal mathematically 
on this basis. 



2.2. Longitudinal Modulation 

The most obvious difference between a true sub-pulse sig- 
nal and a purely periodic signal its longitude-dependent 
amplitude windowing. It is expected that the sub-pulses 
will be scaled in amplitude as a function of pulse phase 
with a shape similar to the average pulse profile. We de- 
note this function a\{(f)), where <j) is the pulse longitude. 

We may also expect some dependence of longitudinal 
sub-pulse spacing on pulse longitude. Expressed another 
way, the sub-pulses of a given pulse are periodic in a non- 
linear function of longitude. The phase must be at least 

^ We denote the apparent period P3 since it may be an 
aliased version of the "true" P3. See Sect. 4.2.1 



monontonic and most likely almost linear, or we would 
not consider the sub-pulses to be (almost) periodic. We 
may therefore treat "sub-pulse phase" [9) in a given pulse 
as the sum of a linear term and a longitude-dependent 
"phase envelope" (pi ((/>)). 

The sub-pulse signal for a given pulse can then be rep- 
resented as 



3(0) = aiA'(0) 



Ak sin[fc(0Pi/P2 +pi((^) + 9') 



(1) 



k=l 



where P2 is an "average" or "nominal" spacing and A' and 
9' are an amplitude scaling and phase shift specific to this 
pulse. The model periodic signal (with a period of P2) is 
composed of a series of Nh harmonics with amplitude and 
phase Ak and 9k for each harmonic k. In most instances 
the sub-pulses will appear almost sinusoidal and the power 
of harmonics with A; > 1 will be small. We therefore ne- 
glect the higher harmonics in further analysis, keeping in 
mind however that in a real signal they may be present, 
with values of P2 and P3 factors of k smaller than those 
of fundamental {k = 1). Through our choice of A' and 9' 
we may define Ai = \ and 0i = so that they may be 
neglected in the case of the fundamental. 

Since the deviation from pure periodicity is repre- 
sented by a longitude-dependent amplitude scaling and 
phase shift, we may simply define a complex modulation 
envelope m\{(j)) = ai((/)) exp(jpi ((/>)), writing the varying 
sub-pulse component of a given pulse as 



s(0) = Rc[mi{(j))A'e"l'P^/P^+^'] 



(2) 



where Re(a;) denotes the real part of x. However, it is 
never the case that sub-pulses subtract from the overall 
emission; if the sub-pulses are to be modelled with sinu- 
soids, we must add a unity bias which itself is amplitude- 
modulated by a\{(j)). In addition there may be some emis- 
sion not associated with sub-pulses (the "non-drifting" 
component, u{(j)))^ which we include here under the simpli- 
fying assumption that it does not vary from pulse to pulse. 
Together these terms are referred to as the "steady" part 
of the emissiorj^. The full signal in a given pulse can be 
written: 

i{<P) = Re [mi(</))A'e**^i/^^+'''l + A' [ai(<^) + u{c^)] . (3) 



2.3. Time Modulation 

In addition to the longitude dependence covered in the 
previous section, drifting sub-pulse signals vary on time- 
scales longer than one pulse period (Pi ) . We define a time- 
dependent amplitude scaling at{t) which gives rise to the 

^ This is not to say that the steady component is constant 
with time, only that it is independent of sub-pulse phase, that 
is its form from pulse to pulse varies only by a pulse-specific 
amplitude scaling that applies to all longitudes. 



R. T. Edwards and B. W. Stappers: Drifting Sub-Pulse Analysis Using the 2DFS 



3 



pulse-specific A' factor used earlier. To describe the time- 
dependence of the spacing of sub-pulses in time, we may 
take a similar approach and define Pt(t) as the difference 
in phase between the true signal and a uniform sinusoid of 
period P3. This function is the source of the pulse-specific 
9' terms of the previous section. Again we can define a 
complex envelope mt{t) = at{t) eKp(ipt{t)) . 

2.4. Full Model of Intensity Signal 

From the preceding two sections we are able to assemble a 
model for the observed intensity of pulsar signals as a func- 
tion of time and pulse longitude. Under the assumption 
that the time-dependent and longitude-dependent modu- 
lations are independent, we may model the overall mod- 
ulation as the product of the derived longitude and time 
windowing: m(</>, i) = mi{(j))mt{t) . Hence the model of the 
sub-pulse component is 



(4) 



where t is taken as constant within a given pulse period. 
We note that in this formulation, if P2 and P3 have the 
same sign, the sub-pulses appear to arrive earlier with each 
successive pulse. We adopt this convention and choose to 
denote the opposite sense of drift with a negative value 
for P3. The full signal including steady components is: 



i{<j>,t) = at{t)[u{(P) + ai{<j>)] 

+ Re [m(0, t)e'(^Pi/P2+2.t/P,)] (5) 

3. The Two-Dimensional Fluctuation Spectrum 

3.1. Definition and Properties 

The two-dimensional Fourier transform of a function 
f{x, y) is defined as follows: 



00 /"OO 



/(a;,2/)e~2"("^+^!^)da;d?/. 



(6) 



— 00 ^ — 00 



where !F denotes the Fourier transform operator. Its in- 
verse differs only in the sign of the exponent: 



basic unit of this decomposition can be seen from the in- 
verse response to a shifted delta function {5(u—uq, v—vq)): 



T-\5{u -uo,v- vo)]{x, y) = e2^^("«^+"o2'). 



(8) 



The form of the real component of this is simply a sinusoid 
in a linear combination of x and y, that is a periodic func- 
tion with maxima and minima lying on lines of constant 
uqx + voy. 

For a discretely, uniformly sampled function with 
Nx X Ny samples, the two-dimensional Discrete Fourier 
Transform (DFT) may be written as 



F{u,v) 



N^-lNy-l 

E E 

j=0 k=0 



fijAx, kAy)e 



— 27ri{uj Ax-\-vkAy) 



(9) 



and is computed over a grid of x Ny values in the range 
-1/2A^ < u < 1/2A^, -l/2Aj, < t> < l/2Ay where A^, 
and Ay are the sampling intervals in x and y respectively. 
As with the one-dimensional DFT, most of the important 
theorems concerning the continuous transform also apply 
to the discrete case. 

When the signal of a pulsar with drifting sub-pulses 
is displayed as a function of pulse number and longitude, 
the patterns made by the sub-pulses and their drifting are 
very similar to those of the inverse Fourier transform of 
a delta function. This fact alone suggests that computing 
the (discrete) Fourier transform of the two dimensional 
signal {i{4>,t)) may be an interesting exercise. 

Since the signal is real- valued, the result of the Fourier 
transform will obey the symmetry I{vi,vt) — 

/(— fj, — where the superscript * denotes the complex 
conjugate. From Eqs. (Q), (H) and (||) the result of this op- 
eration (which we call the Two-Dimensional Fluctuation 
Spectrum or 2DFS) is 



At{iyt)[U{iyi)^ 
+ S{-vi, ~vt) 



AM)] 



(10) 



where functions in upper case denote the (one- or two- 
dimensional) Fourier transform of their lower-case coun- 
terparts. Expanding just one of the two terms arising from 
the sub-pulse modulation s{(f>,t), from the convolution 
theorem we find 



fix,y) = T-'[Fiu,v)] 

F(u, v)e2"("^+^^Mudt>. 



00 pOQ 



(7) 



00 ./ — CKj 



It is usual that the x and y quantities are spatial or 
temporal, whilst u and v are spatial or temporal frequen- 
cies. The effect of the two-dimensional Fourier transform 
is to decompose the input function into a sum of complex 
exponentials, or in the case of a real-valued input function, 
sinusoids of specific frequency, phase and amplitude. The 



Sivi.vt) = T[m\{(l))mt{t)] 

*5[,yi-P,/{2TTP2),ut-l/P3] (11) 
= M{i^i, vt) * 5[,yi - Pi/(27rP2), >^t - I/P3], (12) 

where * denotes the convolution operator. It is also worth 
noting that M{vi, vt) is simply equal to Mi{vi)Mt{vt) since 
the non-zero regions of two responses lie on orthogonal 
lines in the two-dimensional spectrum and the convolution 
reduces to multiplication. 
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So we see that the 2DFS consists of three main compo- 
nents: one for the response of the steady part of the pulsar 
emission, and two mirror-image components for the sub- 
pulse modulation. The steady component is centered at 
DC {vi = i^t = 0), whilst the sub-pulse components are 
centered approximately at vi = fcPi/(27rP2) (cycles per 
radian of pulse longitude), vt — k/P^ (cycles per time 
unit) where k — ±1.^ Since the modulation envelopes are 
not purely real- valued, their Fourier transforms are not ex- 
pected to be entirely symmetrical, and there is no unique 
"center" of the components, however in most cases the 
phase rotation of the envelope will be quite minor and a 
peak around the position specified will be seen. 

The overall width and shape of the components de- 
pends on the nature of the modulation envelopes. For a 
pulsar with no nulling or scintillation, the steady com- 
ponent will have zero width in the I't direction^, and 
a one-dimensional inverse Fourier transform of /(;//, 0) 
in the vi dimension yields the average pulse profile. If 
nulling and/or scintillation are significant, there may be 
some broadening in the vt direction however very frequent 
nulling or scintillation on very short time scales would be 
required to produce a strong, significantly extended com- 
ponent. In addition to the effects of nulling and scintil- 
lation, the sub-pulse components are broadened due to 
phase (or P3) variation over the course of the observation. 

Broadening of the components in the vi axis derives 
from the finite width of the amplitude envelopes (ai((/)) 
and u{(j))) and, in the case of the sub-pulse components, 
due to phase (or P2) variation across the pulse. The width 
of this component is generally quite large due to the small 
duty cycle of pulsars, and will tend to make the spectrum 
significantly non-zero at i/; = 0, = k/P^ if the number 
of sub-pulses observable in each pulse is less than two. 

As noted earlier, the frequency extent of the DFT is 
limited by the inverse of the sampling interval(s). Co- 
efficients for frequencies outside the range may be com- 
puted but will always be identical to those at some (cal- 
culable) point in the range due to aliasing. In the case of 
the 2DFS, the longitude sampling may be chosen arbitrar- 
ily (within the limits of the hardware and the bandwidth 
of the received radiation), but the "time" or "pulse num- 
ber" sampling is naturally set at a value of Pi. For this 
reason, drifting sub-pulse components outside the range 
—0.5 < P1/P3 < 0.5 are indistinguishable from those 
within it. 



3.2. Relationship to the Longitude-Resolved 
Fluctuation Spectrum 

It is readily apparent from Eq. (|^) that the discrete two- 
dimensional Fourier transform can be computed by first 



^ If the sub-pulses are significantly non-sinusoidal, compo- 
nents at other integer values of k will be present. 

* Naturally in practical analysis of data of limited duration 
the time windowing leads to an approximately sinc^ response 
on a frequency scale of 1 bin in the DFT. 



computing a one-dimensional DFT along each column of 
the input data, and then computing the one-dimensional 
DFT along each row of the result (or alternatively, along 
rows first and columns second). The result of the first 
half of this operation when performed on the pulsar sig- 
nal i{(t>,t) is known as the longitude- resolved fluctuation 
spectrum (LRFS; [Backer 1970"a| ,|b|). 

The expected form of this spectrum is easy to under- 
stand in the light of the preceding discussion. Each column 
(i.e. line of constant (f>) in the spectrum will have three 
components, as in the 2DFS. The steady component sim- 
ply appears as a DC term equal to u{4>) + a{4>), whilst 
the two components for the sub-pulse modulation appear 
at zbl/Ps. Calculating the spectrum of both components 
(denoted Ii and 12), we find: 

/i(0, i^t) = mi(0)e^"*^i/^^ X Mt{t) * S{i^t - I/P3) (13) 



and 



. i^t) = "ii((/))e" 



-l0Pl/P2 



X Mtity ^Sii^t + l/Ps) (14) 



That is, the amplitude (or power) spectrum in each phase 
bin is identical excluding the scaling factor ai(0), whilst 
the complex co-efficients are rotated from bin to bin by 
an amount (/1P1/P2 +pi{(f>). 

Due to the presence of a component in both the posi- 
tive and negative half of the spectrum, it is impossible to 
determine the sign of P3 from the LRF power spectrum. 
Signals are effectively aliased to the range P1/P3 S [0, 0.5]. 
Only by examining the sense of phase rotation of the com- 
plex co-efficients as a function of longitude and comparing 
this to the sign in vt of the component being considered 
can the sense of the drift be determined. 

Since the phase relation between longitude bins is ex- 
pected to be quasi-periodic, it makes sense to perform 
a Fourier transform across each row of the complex spec- 
trum in order to determine the sense of the drift and its pe- 
riod (i.e. P2/P1), and to concentrate the sub-pulse power 
in a smaller region of phase space for greater signal-to- 
noise ratio in P3. The result of this operation is the two- 
dimensional DFT of the signal, in other words the 2DFS. 

3.3. Relationship to the Harmonic-Resolved 
Fluctuation Spectrum 

Recently a new technique for sub-pulse analysis has 
been devised by Deshpande fc Rankin (2001 ). It involves 
the computation of the so-called Harmonic-Resolved 
Fluctuation Spectrum (HRFS). Working from the original 
time-series i(t), they compute the one-dimensional spec- 
trum and "stack" it about the v = k /Pi harmonics of the 
signal. Thus, 



H\i]{x,v)^T\i] 



vVPi 



(15) 



where Ti denotes the HRFS and x and y are the fractional 
and integer components of the corresponding spectral fre- 
quency multiplied by Pi . The HRF power spectrum is usu- 
ally plotted with the x parameter as the abscissa (with a 
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domain of to 1) and the y parameter as the ordinate 
(with integer- valued bins up to the extent allowed by the 
sampling interval). 

We note that the 2DFS is in fact intimately related 
(after scaling and rotation) to the HRFS. The former is 
formed by stacking the time-domain data about the pulse 
period Pi and then Fourier transforming, whilst the latter 
is formed by Fourier transforming and stacking the result 
about 1 /Pi . That the two should be essentially identical is 
not immediately obvious but is easy to show (see appendix 

In fact, the decomposition of the one-dimensional DFT 
into a two-dimensional DFT is the basis of parallel FFT 
codes ( pooley fc Tukcy 1965| ). 



We note that values of ly giving rise to 0.5 < Vt < 1 
are aliased to i^t — 1 in the 2DFS (i.e. the opposite drift 
sense), however this is simply a matter of convention: by 
the same notion, signals which truly show the opposite 
drift sense (that is, negative P3P2) are aliased in the HRFS 
into the region 0.5 < ft < 1. In our view, given that there 
is no reason to prefer one sense of intrinsic sub-pulse drift 
over another, specification of P3 in the interval [—0.5,0.5] 
makes more sense than taking a preferred drift direction 
and mapping all drift to this interval (e.g. [0:1]). 

Whilst the two spectra are equivalent, a knowledge of 
their basis in the two-dimensional Fourier transform of the 
longitude-time data is of value in further analysis, as we 
shall see below. 

4. Sub-Pulse Analysis with the Two-Dimensional 
Fluctuation Spectrum 

With an awareness of the mathematical basis of the 2DFS 
in relation to the form of the drifting sub-pulse signal 
described in the preceding sections, a range of analysis 
strategies can be developed for the investigation of partic- 
ular aspects of the sub-pulses. 

The technique was originally conceived with the aim of 
detecting the presence and basic parameters (P2 and P3) 
of sub-pulse drift in pulsars with limited instantaneous 
signal-to-noise ratio. The sensitivity of this method de- 
pends only on the coherence of the drift and the ratio of 
its amplitude to that of any non-drifting signal. The sig- 
nificance of a detection can by improved arbitrarily (for 
stable drifters) by simply integrating more pulses, unlike 
in traditional single-pulse studies where each individual 
sub-pulse must be detectable above the noise. This makes 
it well-suited to studying the sub-pulse properties of a 
large sample of weaker pulsars. 

The complex 2DFS is also of use in studying the details 
of sub-pulse emission in bright pulsars. In this section we 
show how the complex modulation envelopes can be ex- 
tracted from the spectrum, giving full information about 
the deviation from purely periodic sub-pulses. 

4.1. Extraction of the Modulation Envelopes 

In section ^ we saw that the components arising in the 
2DFS due to the sub-pulse drifting take the form of the 



Fourier transform of the combined modulation envelope 
m{(t>,t), shifted to vi = P^'^PiI2tt, vt = Pg^^ To obtain 
the envelope from the spectrum, we begin by masking out 
the steady component^ and the mirror- image component. 
The inverse transform of this spectrum contains a complex 
version of the sub-pulse component of the data which, in 
the limit of perfect masking of other components, is its an- 
alytic signal. Its amplitude is only half of that present in 
the original data due to the loss of the mirror-image com- 
ponent, so to circumvent the effect of this in further anal- 
ysis we multiply the complex co-efficients of the spectrum 
by two after performing the masking. Before performing 
the inverse transform we shift the spectrum by the ap- 
propriate amount horizontally and vertically to place the 
nominal centroid of the sub-pulse response at DC. The in- 
verse Fourier transform of this shifted spectrum gives the 
two-dimensional modulation envelope m{4>,t), as used in 



Sect. 2.4 



In performing the described procedure, is vital that the 
mirror image sub-pulse component and the steady compo- 
nent are both carefully removed from the spectrum. In the 
case where the conjugate mirror image components are not 
separated by regions of nearly zero power, it is likely that 
the Fourier transform of the "true" complex signal had 
power in opposite quadrants, which as a result of the real 
sampling is added to that of its conjugated mirror image 
in the computed spectrum. Under such circumstances, it is 
not possible to reconstruct the true complex signal of the 
sub-pulse component, and the derived envelope will differ 
from that which we had hoped for to a degree dependent 
on how significant the overlap is. The same caveat applies 
if the drifting component is not clearly separated from the 
steady component. 

The two dimensional envelope produced using the 
above procedure can be decomposed into longitude- and 
time-dependent parts under the assumptions made earlier 
about the form of drifting sub-pulses, as such: 



m(0, t) — mi(0)mt(t) + noise. 



(16) 



We use an iterative scheme to perform this decompo- 
sition. By assuming a form for m\{(j)), one can compute 
mt(t) through rearrangement of Eq. (16). Using the re- 
sultant mt(t) SL new mi{(j)) is computed, and so on. The 
process is repeated until convergence, constraining the sys- 
tem to a stable solution by normalising mt{t) to a mean 
amplitude of unity^. Use of simple inversion leads to dom- 



^ This could potentially be achieved by subtraction of the 
mean profile from every pulse, however in many cases the power 
this removes at non-zero frequency (due to the finite time win- 
dowing) is only a small fraction of that present, due to pulse- 
to-pulse intensity variations. 

^ There is also a degree of freedom available for arbitrary 
phase rotation, however the iterative scheme tends to converge 
without the imposition of a constraint to remove this freedom. 
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ination of the inferred phases by noise, so we use instead 
the zero-lag of the normalised cross-correlation function, 



mt(t) = 



Ef:o'm(jA(/),^)mi(jA0)* 



(17) 



(or the equivalent for mi(0)), which is simply a weighted 
version of the longitude-averaged inversion of Eq. ([l6|). 

With noiseless data and adequate removal of unwanted 
spectral components, a convergent result of this algorithm 
clearly gives the correct answer. Under the presence of 
noise, however, care must be taken in the interpretation of 
derived amplitude values. Denoting the RMS noise in the 
real and imaginary parts of the two-dimensional complex 
envelope by an (which may be calculated from the variance 
of the off-pulse samples in the original data) , the variance 
of the noise in the real and imaginary parts of m\ is 



(18) 



(A similar relation applies for a^^: both are reduced from 
fj^ by the integrated power of the complementary enve- 
lope.) 

Taking the amplitude envelope incurs some bias due to 
the presence of this noise. The derived power values are 
over-estimated by, on average, 2af, leading to significant 
bias in parts of the envelope with low (or zero) signal- 
to-noise ratio. This bias is intrinsic to the amplitude- 
phase representation of noisy data, and so is also familiar 
from the estimation of the longit ude-dependent envelope 
from LRF spectra (Backer 1973 ) and of modulation in- 
dices (e.g. Taylor fc Hugucnin 1971 ). A second-order effect 
arises from the fact that the iterative algorithm described 
above includes the noise bias in its unity normalization of 
the time-dependent envelope, leading to less correlation 
and hence an unintended attenuation in the next estimate 
of the longitude envelope. The latter effect is easily ac- 
counted for by correction of the normalization to exclude 
the noise contribution. Compensating for the former effect 
requires a more careful approach. 

At the cost of increasing the variance of the noise, the 
actual (noisy) values of an envelope may used to determine 
the bias to subtract from each element. Alternatively one 
may take — 2a^ as an unbiased estimator of the power 
envelope, however a problem remains when amplitudes are 
to be calculated from negative power values (such as oc- 
curs when the signal-to- noise ratio is low). Nevertheless, 
this is the approach that is usually employed in the calcu- 
lation of modulation indices. A final option that may be 
of practical use is to make a model longitude-dependent 
phase envelope Pm{4') (since in most cases the measured 
phase envelope will be more smoothly varying than the 
amplitude envelope). Neglecting errors in the model en- 
velope, the real part of the product m\((j)) exp{—ipmi(t>)) 
(corresponding in the complex plane to the projection of 



the measured envelope on an axis aligned with the model 
phase for that longitude) should be free of biasQ. 

Once the decomposition is performed, four major 
classes of questions can be asked: 

1. How does sub-pulse phase (or equivalently, P3) vary 
with time? Are there any systematic variations and if 
so what causes them? 

2. How does sub-pulse amplitude vary with time? 

3. How does sub-pulse phase (or equivalently, P2) vary 
with pulse longitude? 

4. How does the sub-pulse amplitude window compare to 
the mean profile? Is there a non-drifting component?^ 

All of the questions listed above can be (and in 
some cases have been) answered to a degree using LRFS, 
TRFS0 or HRFS by taking one-dimensional "stacked" 
power spectra (for amplitude windowing phenomena) or 
one-dimensional "slices" (for examining phase relations). 
However both these techniques lose sensitivity, in the first 
case by summing incoherently and in the second case by 
using less than the full available modulation power. We 
believe that the technique presented here is likely to re- 
sult in measurements of equal or (usually) better quality 
since (as long as the assumptions made are valid for the 
pulsar) the two envelopes form a complete representation 
of the drifting component of the signal. 

4.2. The Carousel Model 

4.2.1. Expected Modulation Envelopes 



The most popular conceptual framework for understand- 
ing drifting sub-pulses is that of a carousel-like system of 
"sparks" uniformly spaced about a ring of constant mag- 
netic elevation and altitude from the star surface, with the 
system as a whole slowly rotating about the magnetic axis 
( Ruderman 1972 ). As the pulsar rotates and its beam pe- 
riodically intersects the observer, a series of pulses, each 
composed of a series of sub-pulses is seen. 

By treating the problem as if all radiation originates 
on the surface of a "polar cap" , and is beamed in a direc- 
tion directly opposite to that of the center of the star (as 
shown in Fig. pi) , we may consider the rotational geometry 



^ This estimator also has the virtue of having Gaussian 
statistics, with half the variance of that which uses the modu- 
lus. 

* If there is no non-drifting component, the longitude- 
dependent amplitude envelope should be identical to the mean 
pulse profile, see Eq. Any difference can be explained by 
the presence of a non-drifting component to the mean profile. 

^ We call the spectrum formed by performing DFTs along the 
pulse longitude axis the Time-Resolved Fluctuation Spectrum 
(TRFS). It could be used to examine sub-pulse phase as a func- 
tion of time, albeit with greatly reduced sensitivity compared 
to the techniques presented here, by examining the phases of 
the coefficients as a function of time for a given frequency (vi) 
bin. 
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to effect, over the course of each successive pulse, a sam- 
pling of the polar cap emissivity along a locus defined by 
the intersection of the sight line and the cap. This simplifi- 
cation gives a means for good conceptual understanding of 
the origin of the observed pulse morphologies and is widely 
used for this purpose (e.g. Lyne & Manchester 1988|; |Mitra 



& Des ipande 1999; Deshpande & Rankin 2001), but it 
must be borne in mind that in reality the radiation is 
expected to be beamed along a tangent to the magnetic 
field line at the emission location. Therefore the observed 
emission at a given instant relates to a plasma fiux tube 
with its foot at some point on an arc on the polar cap con- 
necting the depicted "line of sight" and the magnetic pole 
(labelled p in Fig |l|c). However, since the field line tangent 
angle at a given altitude (where the emission originates) is 
approximately proportional to its magnetic colatitude an- 
gle, the derived sight-line locus differs from the true polar 
cap sampling by a simple scaling of the magnetic colati- 
tude co-ordinate (see also Lyne k, Manchester 198^ )^ In 
any case, the true sampling effected in magnetic azimuth 
is the same as that derived assuming radial beaming, and 
if the polar cap emission pattern is a ring of sparks (or ra- 
dial "spokes" ) , it is this co-ordinate alone that determines 
the pulse longitude of the observed sub-pulses. 

In systems where the sight-line makes a tangential pass 
of the ring, the emission from the sparks is seen as a se- 
quence of one to a few almost equally-spaced sub-pulses, 
the phase of which "drifts" in time with respect to a fidu- 
cial point in pulse longitude due to the rotation of the 
carousel. The number of sparks present, and the viewing 
geometry and the angle between the spin and magnetic 
axes determine the longitudinal spacing of the sub-pulses 
(P2), whilst the time spacing (P3) depends on the rotation 
rate of the carousel and the number of sparks it contains. 

There are two main effects predicted under the 
carousel model that are of importance here. 

Firstly, the carousel rotates continuously (with some 
period P4 = where there are N sparks), meaning 

that any variations in spark intensity, shape or spacing 
about the ring should manifest as a periodic feature in the 
drifting sub-pulse signal. This periodicity would appear as 
a periodicity in the time-dependent amplitude and phase 
envelopes, giving rise in the 2DFS to a pair of "sidebands" 
at frequencies I/P4 higher and lower than the primary 
feature at i^t = I/P3. This was the startling finding of 



Deshpande fc Rankin (2001) who applied the HRFS for 



The second effect predicted under the carousel model 
is that the longitudinal spacing of the sub-pulses depends 
on pulse longitude. The sparks are presumed to be spaced 
uniformly in magnetic azimuth [tp) and their separation in 
spin longitude is known, but the mapping from pulse lon- 
gitude to magnetic azimuth for the sight line is dependent 
on the degree of spin-magnetic pole misalignment (a) and 
the viewing geometry (where C is the angle made between 
the line of sight and the positive spin pole). The relation 



tamp 



sin (/) sin ^ 



cos C sin a — cos </> sin ( cos a ' 



(19) 



the fir^ t time, to observations Of PSR B0943+10. The dc- 



where the signs of the numerator and denominator deter- 
mine the quadrant for ^p in the usual way (see appendix 
^ for derivation) , and (f> is measured relative to the pulse 
longitude at which the sight-line makes its nearest traverse 
of the sight-line. 

The mapping of pulse longitude to magnetic azimuth 
is monotonic in the region in which emission is seen, with 
a gradient of the opposite sign to sgn/? = (where 
/3 = ( ~ a). In our analysis, P2 and hence the rate of 
change of sub-pulse phase {d9/d(p) are always positive, so 
a factor of — sgn/? will appear in the expression relating 
0{(f)) to magnetic azimuth (ip((p)). 

The fact that the carousel itself rotates with respect 
to a fixed point of magnetic azimuth must also be taken 
into account. The rotation period of the carousel is 27r/a;p 
(see fig. 1^) where positive ujp indicates counter-clockwise 
rotation viewed from above the active magnetic pole. The 
direction of drift seen in the sub-pulses depends also on 
sgn P : if their signs are opposite (as is the case in fig. |lj) 
the sub-pulses will drift toward the trailing edge of the 
profile. We label this with a negative value for P3, via 
P3 =27r(sgn/3)/(7Vt^p). 

Since the sub-pulse phase at a given longitude is only 
sampled once per rotation period (Pi ) , we will always ob- 
serve some periodicity —0.5 < P1/P3 < 0.5. Under the 
carousel model it is possible that each spark drifts in longi- 
tude by more than half of P2 from one pulse to the next. In 
this case we may say that the observed P3 is an "aliased" 
version of the true time taken for the carousel to drift by 
1/iV turns, P3. The ambiguity in the "true" P3 can be 
expressed in terms of the (signed, nearest whole) number 
of sparks (n) that drift by undetected from one pulse to 
the next: P1/P3 = n + P1/P3. As shown by Deshpande 



tectability of tlii ij effe ct UHing 2DFS Leclmiquetj in ditjcutjjjcd 
further in Sect. 4.3.3. 



Furthermore, we note that re fraction in the magnetosphere 

(e-l 



Petrova & Lyubarskii 2000 ) would also alter the sampling 



& Rankin (2001), this aliasing can potentially be con- 
strained with a measurement of P4 by observing that 
P4/P1 = \N/{n + Pi/P3)\ and solving for integer pairs of n 
and N (in their case 37.35±0.52 = 7V/(?i-0.4645±0.0003) 
implying n = 1 and N = 20 f^ 



of the polar cap afforded by the sight line. Again, only the mag- 
netic colatitude of the derived emission location is corrupted 
(assuming a spherically symmetric plasma distribution), but 
in this case the mapping from ray angle to emission location 
is not a simple scaling and in some cases may not even be 
one-to-one. 



Deshpande fc Rankin (200l[ ) cited the observed longitudi- 



nal dependence of the p olariz ation position angle as support for 
thi s solution (see Sect. 4.2.2). A lthough not specifically noted 
by Deshpande & Rankin (2001), the polarization constraint is 



in fact vital to the argument: without it other aliasing solutions 
(all with A'' > 54 and |n| > 1) are permitted. 
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The expression for 9{(j)) is the sum of three terms, the 
first a sign-corrected and scaled version of ip{<p) (Eq. |l9| ), 
the second a term which accounts for carousel rotation 
over 
sentine 



th|e course of the pulse, and the third an offset repre- 
the carousel position at 



independent of the instantaneous P3 ( Biggs et al. 1985 ). 
Other pulsars that show major changes (but not sign re- 
versa ls) in the drift rate are PSR B201 6+28 ([Taylor et al. 
1975|) and PSR B0809+74 (after nulls; [Page 1973| ; |Lyne fc 



0. 



6l(</)) = -iVsgn/3tan"^ 



sin (f> sin C 



cos C sin a — cos (p sin C cos a 



Ashworth 1983), both of which also show constant longi- 
tudinal spacing. It therefore appears safe to assume that 
most systems are not highly aliased and the longitudinal 
spacing can be assumed constant with time. 

(20) 4.2.2. Comparison with Polarization Behaviour 



Near the point of zero longitude, the rate of traverse 
of magnetic azimuth is at its maximum: 



(21) 



The apparent rate of traverse of sub-pulse phase is 
therefore given by 



dV- 




sinC 


d0 


max 


sin/3 





sinC 






sin f3 



Pi 

A' 



(22) 



and will generally be close to the "nominal" frequency 
used in 2DFS analysis (via AO/A^ = PilP^)- The differ- 
ence between the predicted phase evolution and a con- 
stant slope of this gradient gives the predicted longitude- 
dependent phase envelope. An example is shown in fig. 
^. We note that the longitude-dependent phase envelope 
arising from a rotating carousel will always have a similar 
form (given sufficient signal to noise ratio). The gradient 
of 0(0) is at a maximum at (/) = 0, and declines symmetri- 
cally on either side. The spectral shift, if the correct value 
of P2 is chosen (and this may be fine-tuned by examining 
the resultant envelope), will result in an envelope with zero 
gradient at </) = 0, and negative gradients at all other lon- 
gitudes, with the steepness increasing towards the leading 
and trailing edges of the profile. 

The techniques used in this work decompose the sub- 
pulse signal into the product of two complex envelopes 
which depend only on time and longitude respectively. 
As such, it is important that the longitudinal spacing 
of the sub-pulses does not vary with time, and that the 
time-spacing of sub-pulses does not vary with longitude. 
The longitude-independence of P3 is trivially true in the 
carousel model, since the carousel rotates as a whole and 
the angular speed is constant around the ring. This is seen 
in LRF spectra of virtually all pulsars with drifting sub- 
pulses, see e.g. Backer (1975 ). The time- independence of 
P2 is technically not satisfied under the carousel model, 
since the longitudinal spacing depends on P3 (see Eq. ^ , 
which in turn depends on the circulation rate, which may 
vary with time. However, for this effect to be significant 
Eq. ( p^ ) must be dominated by the carousel rotation, i.e. 
|n| ^ 1. In this case, even minor variations in the rotation 
rate could cause the observed (highly aliased) P3 to vary 
wildly in sign and magnitude. To our knowledge, the only 
pulsar to clearly show sense reversals in its drift is PSR 
B0826— 34, and the longitudinal spacing is (within errors) 



The longitudinal dependence of the observed position an- 
gle of linearly polarized radiation is expected to be inti- 
mately related to carousel rotation sub-pulse behaviour, 
and comparisons between polarimetric and drifting-sub 
pulse analyses have the potential to provide more infor- 
mation than either alone. Under the "magnetic pole" or 
rotating vector model" (RVM; Radhakrishnan & Cooke 
1969), the position angle is that of the sky projection of 



the magnetic field line on which the emission region at 
a given instant lies. Analogous spherical geometry to that 
used in the derivation of the magnetic azimuth of the sight 
line leads to a similar relation: 



tan —X ■ 



sm (p sm OL 



cos a sin C — cos </> sin a cos C 



(23) 



(See appendix |B| . Eq. (^3|) differs to the formula of 
pComesaroff 197C and virtually all later works on pulsar 
RVM fitting in that it assumes the standard astronomical 
convention of position angle increasing counter-clockwise.) 

In principle, a and ^ can be determined directly by 
fitting the observed position angle swing. These values 
could then be used with a measurement of the longitude- 
dependent sub-pulse phase envelope and Eq. (|2^) to de- 
termine iV and check the validity of the carousel model 
(both through its ability to fit the data and by unambigu- 
ous consistency of the fitted value of N with an integer). 
In practice there may often be little measurable curvature 
in either function, in which case it makes more sense to 
work with the (fitted and presumed constant) gradients of 
the curves. 

The first thing to notice is that the position angle 
swing may have either positive or negative slope, depend- 
ing on the sign of /3. The maximum rate of position an- 
gle swing occurs at the closest point to the magnetic 
pole, where dx/d0 = — sin a/ sin /3. If the degree of sub- 
pulse aliasing is known (or presumed zero) we can de- 
termine the physical direction of carousel circulation via 
c^p = 27r(sgn/3)/(A^P3). 

Secondly, the difference in the absolute value of the 
slopes can provide useful information. Since sin ('/sin/? — 
sin a/ sin /3 ~ cos(q; -f- /3/2) (for small /3), one may take 



(24) 
.(25) 
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Unless the degree of aliasing is very high (with of the order 
of one complete carousel rotation for every star rotation) , 
the second term may be neglected since it is much smaller 
in magnitude than the likely uncertainty in cos(a + /3/2). 
The first term may be estimated for example from the 



pulse width (and position angle slope; Lyne & Manchester 
1988|~plankin 1990), or taken as unknown to yield N 
(d0/d(^)/(|dx/d0|±l). In certain fortuitous circumstances 
use of this expression may determine N to better than 
±1, and we could use the nearest integer to calculate a 
and (3 using the known position angle and sub-pulse phase 
slopes. 



We note that the work of Dcshpande fc Rankin (2001 ) 



arrived at d^/d0 — dx/d(j) — — sgn/3[U, seemingly under 
the assumption of a ~ 180°. Their result of iV = 20 for 
PSR B0943-f 10 remains valid under the analysis of the 
preceding paragraph. Taking their value of 34 for the max- 
imum sub-pulse phase slope, using the position angle slope 
of —2.5 from ^uleymanova et al. (198g ) and neglecting the 
(probably significant) uncertainty in this measurement, we 
find N — 14 ± 6. The marginal consistency may indicate 
that cos(q; -I- (3/2) is close to —1, implying that a is close 
to 180° and the neutron star rotates clockwise as viewed 



by the observer, in agreement with Dcshpande & Rankin 
(200l|)] If measurements, with quoted uncertainties, of the 
position angle and sub-pulse phase slopes with quoted un- 
certainties were available, it would be possible to use the 
TV = 20 result to constrain a and (3 without the usual 
recourse to the use of pulse width-period relations. 



4.3. Tests on Simulated Data 



Whilst the techniques outlined in Sect. [4.l| are simple and 
quite direct, we feel that it is important to demonstrate 
that they perform as expected. To this end, we used the 
carousel model to produce several simulated time series 
which we subjected to 2DFS modulation window estima- 
tion. 



4.3.1. Form of the Model 

For each sample, we computed the pulse longitude (0), and 
the corresponding magnetic azimuth (ip) that is sampled 
by the sight-line at that longitude using Eq. (|l^) and a 
geometry specified by a and f3. This was converted to a 
"sub-pulse phase" (9) through multiplication by — sgn j3N, 
and addition of an offset 9d{t + (j>Pi/2Tr) (recalling that t is 
taken as constant in any given pulse), advancing at a rate 
of one turn per interval P3 (t) , which produces the drifting 
as a function of time. To investigate the detectability of 
variations in P3 , we simulated the effect of a pulse null as 
a period of no drifting (P3 — > 00) during the null, followed 



by an exponential recovery to the asymptotic rate P3 . The 
full form of Psit) is given by the following equation: 



P' 



P^l-e-*/^)- 



t < tn 
t>tn+d 



(26) 



where tn is the time at which the null begins, d is its 
duration, and r is the time constant for drift rate recovery. 

For each sample the flux contribution from each spark 
was computed as a Gaussian in magnetic azimuth (of 
width <Js and a peak value of unity) and the contributions 
added. The result was scaled by a factor of 1 + csmij^ 
(0 < c < 1) to simulate an asymmetry in the azimuthal 
flux distribution of the spark carousel, to test for de- 
tectability of the carousel circulation time. 

Longitudinal amplitude windowing was applied 
through multiplication by a Gaussian in cj) (of width am, 
centered at = with unity peak), and the resulting 
intensity was in turn scaled by a factor 



1 t<tn 
at{t) = { tn<t<t„ + d 

1 t> tn+d 



(27) 



to simulate the switching-off of the emission during the 
null. To this was added a random noise value (drawn from 
a Gaussian distribution of zero mean and standard devia- 
tion cr„), giving the value i{4>,t) to use in later analysis. 

This model leads to a sequence of pulses, each com- 
posed of sub-pulses whose amplitudes vary (in the spe- 
cific simulations below) from ~ 10% to ^ 100% of the 
unity-peaked overall longitude envelope. The modulation 
envelopes measured from the simulated results are com- 
pared with a predicted envelope with amplitudes given 
by half the peak-to-peak sub-pulse variation produced by 
the noiseless model at a given longitude, after normaliz- 
ing for a unity-mean time-dependent amplitude envelope. 
The predicted longitude envelope peaks at an amplitude 
of ~ 0.45, compared to ~ 0.50 for the average profile: 
that is, a small non-drifting component is present, since 
the model carousel has some emission at all points on the 
ring. 

4.3.2. PSR B0809-F74-like Model 

We chose a set of model parameters that give a pulse se- 
quence similar to those observed from PSR B0809-I-74 at 
frequencies around 400 MHz. Specifically, w e used a ge- 
ometry of a = 9°, /3 = 4.5° ( [Rankin 199"3a| |b|), requiring 
(from Eq. ^) = 16 to reproduce the reported P2/P1 
of 40 ms/1.29 s ~ 32. A value of U.OPi was used for P^. 
The main pulse window had a full width at half-maximum 
(FWHM) of = 12° (fromcr™ = 0.014 turns), with each 
sub-pulse having a FWHM of 5.0° (from = 0.2/ N). 
With the exception of a and /3, all parameters above 
derive from Taylor ct al. (1975). We included a null at 



Deshpande and Rankin used a different convention, where 
a 180° - a and /3 ^ -/3. 



tn = lOOPi, of duration d = 8P1 and recovery time con- 
stant T = 12Pi. No spark-to-spark intensity variation was 
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included in this model (i.e. c = 0) since a time-varying 
drift-rate complicates the analysis of this behaviour]^. We 
produced two data sets, one with no noise ((T„ = 0) and 
one with moderate noise levels ((t„ = 0.5), each with 512 
pulses in 4096 longitude bins, of which only the central 
512 were used. 

The noise-free data and the resultant 2DF power spec- 
trum are shown in Figs. ^ and ^ respectively. The form of 
each is as expected from the considerations of the pre- 
ceding sections. In particular, the spectrum consists of a 
DC component, and the first and second harmonics of 
the drifting response, with all components significantly 
extended in the horizontal axis (due to the longitudinal 
windowing) and also at a low level in the vertical axis 
(due to the nuUing-induced time windowing). For clearer 
presentation, the spectrum in fig. ^ was produced with 
enhanced resolution by padding the input data with ze- 
roes before performing the DFT. In practice this would 
usually not be done, as it is unnecessary for the determi- 
nation of a nominal P2 and P3, and complicates evaluation 
of the significance of concentrations of power. We there- 
fore take the position of the coordinates of the peak bin 
of the 2DFS computed at the raw (non-padded) resolu- 
tion as our nominal periodicity, giving P2 — Pi/32 and 
P3 = Pi/0.092. This is in agreement with the parameters 
used to generate the data, given the bin size of 8 x 0.002 
in P1/P2 X P1/P3. 

After forming the full (symmetrical) complex spectrum 
and shifting it by the nominal values for I/P2 and I/P3, 
the responses of the average profile and the mirror-image 
of the fundamental were removed with a notch filter. To 
reduce the degree of "ringing" in the impulse response of 
the resultant spectrum, the filters incorporated a gradual 
transition of the form 

f In-i^tcl/w <Q.5 

fiiyt) = <^ 2\iyt - tytcl/w - 1 0.5 < \iyt - v^cMw < 1 , (28) 

where fiyt) is the filter transmission as a function of fre- 
quency in the Vt axis, vtc is the centre frequency of the 
component to be removed and w is half the total width 
of the filter, in this case w = 0.01/Pi. Since the com- 
ponents to be removed were broad in the vi axis (and 
certainly overlapping in this parameter with the desired 
component, now shifted to DC), the filter transmission 
was independent of frequency in this axis. 

The resultant spectrum was then scaled by a factor 
of two to account for the use of only one of the two 
(conjugate-pair) components over which the modulation 
power is spread, and the inverse Fourier transform was 
taken. Using the scheme described in Sect. 4.1, we then 
decomposed the resultant two-dimensional complex mod- 
ulation envelope into the product of longitude and time 
envelopes. The results are shown in Figs. |^ and ^. 



It is evident from the figures that some systematics 
are present in the residuals of envelopes of both axes. 
Beginning with the time-dependent envelope of the noise- 
less data, we see that there is some low-level periodic vari- 
ation in both the amplitude and phase terms (panels b and 
e), with a period equal to P3. We believe that this is due 
to a combination of two effects. The first is the presence of 
the second harmonic of the sub-pulse modulation, which 
is shifted in the analysis to the former position of the fun- 
damental. The second is the impulse response of the notch 
filters used to removed the DC component. This has a sim- 
ilar form to a sine function (with a somewhat shorter time 
response due to the choice of transfer function), producing 
the enhanced oscillations around the null and the bound- 
aries of the data set. The effect of the latter can be reduced 
by making a smoother transition in the filter, however this 
necessitates the loss of a greater range of the spectrum, 
which has the effect of corrupting both the amplitude and 
phase response to the post-null frequency recovery. In our 
view the phase evolution is more likely to be of interest 
than highly accurate amplitude evolution information, so 
a narrow filter was chosen. 

As can be seen in panel (e), there is very little error 
in the estimated phase, beyond the low-level zero-mean 
P3 periodicity. The residuals of the noisy data appear, 
as expected, consistent with the addition of uncorrelated 
zero-mean noise to the values estimated from the noise- 
less data. The degree of scatter seen here indicates that 
estimation noise is likely to dominate over the systematic 
effects described above for data of all but the highest of 
signal-to-noise ratios. An exception is the time during the 
null itself, where the phase is of course not measurable, 
and the estimated amplitude is affected by the filter re- 
sponse. This error is trivial to correct by examining the 
pulse-to-pulse mean flux density to identify nulls. 

The longitude-dependent modulation window (Fig. |6|) 
also shows some systematic error in the estimates made 
from noiseless data. We believe this is largely due to errors 
in the estimated time-dependent phase envelope discussed 
earlier. This is a second-order effect, and its magnitude 
is very small, in the amplitude window appearing at the 
one percent level. The systematic trend in the amplitude 
window estimated from the noisy data is a simple arti- 
fact of the representation of noisy complex data as ampli- 
tude and phase values, as discussed in Sect. 4.1. Since the 



variance of the noise in the synthesized data was 0.25, 
and given the doubling of spectral coefficients to com- 
pensate for filtering out the mirror image component, the 
mean power of the noise in the 2D envelope is expected 
to be 1.0 per complex coefficient, giving ct^ = 0.5. This 
leads to a variance of ~ 1.0 x 10"'^ in the real and imag- 
inary parts of the longitude-dependent envelope (since 
Sj!!o ^ l™t(jPi)P — Nt = 512), and hence a bias of up 
to ~ 0.045 at the edges of the envelope. Techniques for 



removing this bias are discussed in Sect. 4.1. The second 



We note, however, that after determining the drift rate as 
a function of time, the data could in principle be re-binned 



or interpolated to remove the effect of the variation, leaving a 
pure periodicity from the circulation. 
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order effect of underestimation by a normalization error 
(Sect. O) is expected to amount to only ~ 1%, making 
it insignificant and not visible in the plots. 



4.3.3. PSR B0943+10-like Model 

As an alternative that illustrates the applicability of these 
techniques to other observed phenomena, we produced 
simulated data from a model with similar parameters to 
PSR B0943+10. 

We assumed a geometry of a = 11. 5°, /3 — 5.4° 
([Rankin 19934 H) with = 20 sparks ( |Deshpandc~fc 
Rankin 200 1|). A profile width of 20.5° FWHM ( |Go^ 



fc Lyne 1998| ) (giving Um — 0.023 turns) was used, with 
sparks evenly spaced with widths given by a a = Q.2/N, 
as m Sect. [1.3.2|. A value of Pi /0. 5355 was used for Pg 
(Deshpande & Rankin 2001), with a moderately strong 
dependence on spark number (c = 0.2) applied to give 
rise to a measurable signal of carousel circulation. No nulls 
were included in the simulation, but the instantaneous P3 
values were made to have some variation in order to model 
the observed (large but probably finite) Q (= vjl^v). 
To this end we produced samples from an uncorrelated 
Gaussian noise process and filtered them in the frequency 
domain with the function g(y) = exp[— i^/O.OlHz], to pro- 
duce a slowly varying function which was scaled to have a 
range of ±0.05/P3. This was added to to the no minal I/P3 
to give an instantaneous 1/P3(i). As in Sect. 4.3.2 , one 
noiseless and one noisy ((7„ = 0.5) data set was produced. 
Each data set contained 512 pulses in 2048 longitude bins, 
of which only the central 512 were used. 

The noiseless data and the resulting 2DF power spec- 
trum are shown in Figs. |^ and ^. Notice that the I/P3 
feature of the fundamental (at 0.5355Pi) is aliased to 
a frequency of — 0.4645Pi. This reflects the fact that if 
each sub-pulse is associated with its closest neighbour 
in pulse longitude from the following pulse, its drift ap- 
pears to progress from earlier to later longitudes. This 
approach treats all measured P1/P3 as coming from the 



range ±0.5. The conclusion from the work of Deshpande 



fe Rankin (2001) is that each sub-pulse actually drifts 



0.5355 times the sub-pulse spacing towards earlier lon- 
gitudes with each successive pulse. Hence, the "nearest 
sub-pulse" (IP1/P3I < 0.5) convent ion m akes a sub-pulse 
counting error of n = -|-1 (see Sect. 4.2.1 ). It is important 
to note that this is simply a choice of convention. The 
two cases are indistinguishable without additional clues 
as employed by Deshpande & Rankin (2001). 



The 2DFS were filtered with a function of the form 
shown in Eq. (p8|), with w = 0.03/Pi, and shifted to place 
the observed peak of the fundamental (Pi/P3 — —0.465, 
P1/P2 = 40) at DC. This was scaled by a factor of two 
and inverse transformed and the result decomposed into 
time- and longitude-dependent modulation envelopes as in 
Sect. 4.1 and Sect. 4.3.2, The results are shown in Figs. ^ 
andlloT 



As with the 0809-1- 74-like data, some systematics are 
present in the results from the noiseless data. In this case 
they are less severe due to the good separation between 
the fundamental of the drifting component and the DC 
component. In any real (noisy) data the errors are likely to 
be dominated by noise. In this model, due to the broader 
simulated main pulse, the phase rate (i.e. P2) variation 
across the profile is noticeable even in the noisy data. As 
expected, the bias in the wings of the longitude-dependent 
amplitude envelope is of the same magnitude as that in 
the 0809-h 74-like data. 

The amplitude variations due to the simulated carousel 
circulation are clearly visible in the modulation envelope 
estimated from the noisy data, as are the phase variations 
due to the simulated frequency noise. Since the phase vari- 
ations are quite small, the amplitude modulation appears 
highly periodic and can be detected in the power spec- 
trum of the inferred amplitude envelope. (Fig. |l^). The 
modulation also appears as a pair of "sidebands" around 
the steady and drifting com ponents in the 2DFS (Fig, p ). 
After previous authors (e.g. Deshpande fc Rankin 2001), 



we have "stacked" the columns of the power spectrum 
corresponding to I/P2 < 64 to produce an "average" 
fiuctuation power spectrum. Fig. |l^ shows the stacked 
power spectrum around the — 0.4645P]^^ drifting compo- 
nent. The sidebands associated with the carousel rotation 
are clearly present, offset by 0.5355Pf V20 0.0268Pf ^ 
from their respective parent features. 

The significance of detections of the carousel circu- 
lation measured via the power spectrum of the inferred 
amplitude envelope versus the stacked 2DF power spec- 
trum must be evaluated carefully. Since the former is esti- 
mated from a single Fourier transform, the baseline noise 
power (suitably scaled) is distributed with two degrees 
of freedom. In contrast, the stacked spectrum is formed 
by adding 16 columns of the two-dimensional spectrum, 
so the (suitably scaled) spectral values follow a dis- 
tribution with 32 degrees of freedom. For a detection of 
97.8% confidence (corresponding to the 2-sigma point if 
the noise was Gaussian), the threshold normalised spec- 
tral power values are ^ 3.8 for the power spectrum of 
the inferred amplitude envelope, versus ~ 1.6 for the nor- 
malized stacked spectrum. The corresponding "4-sigma" 
points are ~ 10.4 and ^ 2.3. Neglecting effects due to 
phase noise, each responds to a modulation of amplitude 
c with a spectral power value proportional to c^, plus the 
"noise floor" value of (on average) unity. From the values 
of peaks observed in these simulated data we can infer 
a proportionality constant of ^ 20 between the (noise- 
subtracted) responses of the two types of spectra. Hence, 
values of c giving rise to "2-" or "4-sigma" detections in 
the stacked spectrum (with power levels of 1.6 and 2.3) 
would produce far more significant detections (with power 
levels of ~ 13 and ^ 27) in the power spectrum of the 
inferred amplitude envelope. This enhancement in sensi- 
tivity (amounting to a factor of ^ 2 in the minimum de- 
tectable amplitude of carousel circulation modulation) can 
be understood as a result of exploiting the inherent phase 
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relations of the complex spectrum, rather than (incoher- 
ently) summing spectral power values. 

As noted earlier, the model assumed here also pro- 
duces sidebands around the DC component of the spec- 
trum. Such sidebands were not observed in the studies 
of Dcshpande fc Rankin (200l| ), a fact that needs to be 
understood if carousel circulation is to explain the other 
sidebands. One possible explanation is that the total (inte- 
grated) luminosity of each spark of the carousel is approx- 
imately equal (leading to no sidebands around DC) but 
that the width of sparks differs as a function of position 
in the carousel. 

5. Conclusions 

We have shown that the two-dimensional Fourier spec- 
trum of pulsar longitude-time data is of value in the anal- 
ysis of drifting sub-pulses. We consider the drifting compo- 
nent of the signal as the convolution of a complex "modu- 
lation envelope" with a pure two-dimensional periodicity. 
The modulation envelope describes the variations in aver- 
age sub-pulse amplitude and phase as a function of longi- 
tude and time, and can be decomposed into the product 
of two one-dimensional envelopes which are functions of 
time and longitude respectively. This makes the technique 
well-suited to studying a variety of phenomena associated 
with drifting sub-pulses, including longitudinal variations 
in sub-pulse spacing (induced due to viewing geometry or 
other factors), variations in the drift rate as a function 
of time (due to the recovery from a null, random phase 
noise, etc.), comparison of the longitudinal amplitude de- 
pendence of the drifting and steady components of the 
pulsar emission, and variations in the average sub-pulse 
amplitude as a function of time (due to nulling, carousel 
rotation, etc.). 
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Appendix A: Equivalence of the 2DFS and HRFS 

As noted earlier, the 2DFS and HRFS are in fact equivalent. 
In this section we briefly show how this fact arises. 

Consider the response to a single sinusoid, i{t) = sm2-Kvt. 
This will appear as a delta function in the HRFS at x = 
Frac[!/Pi], y = Int[i'Pi] (both in units of cycles per Pi-interval), 
where Frac and Int denote the fractional and integer parts of 
their arguments. For the 2DFS, when the data are stacked 
the signal will appear in each pulse period (of length Pi) 
as a sinusoid with frequency vP\/2-k (cycles per radian of 
pulse longitude) and phase 27rt/Pi x Frac[i/Pi] where t/Pi is 
the integer pulse number. Transforming first across rows re- 
sults in a delta function at ui = yP\/2-K with a phase an- 
gle of 27ri/Pi X Frac Pi] (and a corresponding component at 
—vP\/2n which can be ignored for now as it is simply the 
source of the symmetry of the 2DFS). Considering the column 



corresponding to = uPi/2n, the response is a pure complex 
exponential with a frequency of Frac[:/Pi]/Pi (cycles per time 
unit). Therefore the result of the two-dimensional DFT is a 
delta function at ui = vP\/2-n, vt — Frac[!/Pi]/Pi. A 1:1 map- 
ping of Ut = a; Pi, — y/2n thus applies in this case, when the 
fact that the fi bin size in the discrete 2DFS is l/27r is consid- 
ered. Since both transformations are linear and 1:1 invertible, 
we therefore conclude that they are identical in all cases, given 
appropriate mapping between parameters. 

Appendix B: Derivation of Magnetic Azimuth 

The derivation of equation |l^ follows a straightforward coor- 
dinate frame transformation, from the spin frame (in which 
the sight line parameters —0 and ( represent longitude and 
co-latitude), to the magnetic frame. In each frame longitude is 
defined such that the positive spin pole of the complementary 
frame is at longitude zero. We definine basis vectors [x y z] 
and [ x' y' z' ] for the spin and magnetic frames respectively, 
aligning z (and z') with the positive poles and placing x and 
x' at longitude zero in their respective frames. In what fol- 
lows all vectors (including d and /i) are taken as being of unit 
length. The expression for d in the spin frame is: 



(B.l) 



The transformation between frames may be accomplished 
by rotating about y by the angle —a (to align the pole with 
the magnetic pole), then rotating in longitude (i.e. about z') 
by TT radians (to place the spin pole at longitude zero). The 
transformation is expressed as follows: 
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Deriving the longitude of d in this frame (i.e. the magnetic 
azimuth) via tan?/) = dyi/d^i yields equation [l£|. This result, 
after differentiation and correction fo r a difference in the sign 
convention for ^, agrees with that of Ashworth (1988). 

The derivation of the polarization position angle follows a 
similar path. It is based in a third frame which is fixed to local 
rest frame the observer, has its origin at the star center, pos- 
itive pole pointing to the viewer (i.e. z" = d) and places the 
positive spin pole at longitude zero. All meridians in this frame 
project as lines on the sky with position angles offset from the 
position angle of the spin axis by their longitude (x). Under 
the RVM, the position angle of linear polarization is given by 
the projection of the local magnetic field line on the sky, which 
is given by the projection of the great circle co-planar with the 
fi and d. Hence the (counter-clockwise) polarization position 
angle (minus the position angle of the projected spin axis) is 
given by the longitude of the magnetic pole in this frame. The 
transformation from the spin frame is a rotation of <j) in lon- 
gitude, followed by a rotation about y by a (versus ^ in the 
above) and a rotation of tt radians about y" as above, whilst 
the input vector has spherical coordinates (0, a) (versus {—(f), C) 
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in the above). From this it follows that equations ^ and 
differ only through switching of a and C, and a sign reversal (in 
or equivalently, in the final longitude value). 
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Fig. 1. Diagram of sub-pulse emission geometry. In part 
(a) the angular momentum, magnetic moment, and line- 
of-sight vectors are shown with symbols h. fi, and d re- 
spectively. The pulsar rotates about the spin axis with an 
angular velocity uji = 27r/Pi, whilst the ring of sparks 
that gives rise to the sub-pulses rotates about the mag- 
netic axis with an angular velocity ujp — —2'n/NP^, where 
= 10 is the number of sparks. The sight-line rotates 
relative to the pulsar spin frame and is shown here during 
the "on-pulse" . Its intersection with the star surface or 
the emitting surface is shown in bold. The co-rotation of 
the magnetic field (and the emission associated with the 
pulse) and the star causes the sparks to appear as sub- 
pulses that are amplitude-modulated at the main pulse 
(i.e. spin) period. Within this modulation window the lon- 
gitudes of sub-pulses drift from one pulse to the next due 
to their physical motion about the magnetic pole. Part (b) 
shows the main vectors as they appear in the plane they 
share when the sight-line makes its closest approach to the 
magnetic pole. The angles C (between the sight-line and 
the spin axis), a (between the magnetic and spin axes) 
and /3 = C — a (between the sight-line and the magnetic 
axis at their closest approach) are shown. In part (a) a 
second line of sight {d') is drawn, along with the merid- 
ian it shares with with the magnetic pole. The spherical 
triangle formed by i, and d' is duplicated in part (c). 
In addition to a and C, this diagram illustrates the mean- 
ing of pulse longitude ((/>), polarization position angle (x, 
measured relative to the position angle of the spin axis) 
and magnetic azimuth [tp) and co-latitude (p). 
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-5 5 

Pulse longitude. (deg) 

Fig. 2. Illustration of the variation in sub-pulse phase (6) 
as a function of pulse longitude (0). We used a = 10°, 
P — 2° and = 10 in this example. The dashed line 
shows the value expected through the spherical geome- 
try (Eq. 20), the dotted line shows that of a pure peri- 
odicity with equal phase slope at = to the dashed 
line (Eq. and the solid line shows their difference. 
This latter quantity can be considered to be a longitude- 
dependent phase-modulation envelope to be applied to the 
constant periodicity to produce the real signal. 
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Fig. 3. Simulated PSR B0809-h74-like data. A sequence 
of 100 pulses from the noiseless data are shown, with the 
interval chosen to show the null. 
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Fig. 4. Two-dimensional fluctuation power spectrum of 
simulated PSR B0809+74-like data. To enhance low-level 
structure, grey-levels have been over-saturated such that 
pure black corresponds to 10~^ times the peak power. The 
resolution used here is higher than that of the direct DFT 
and is achieved by zero-padding the input data. 
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Fig. 5. Complex time-dependent modulation envelope 

and residuals for PSR B0809-|-74-hkc simulated data. The 
central panels (c and d) show the amplitude and phase re- 
sponse expected from the simulation parameters and the 
chosen nominal P2 and P3. Panels (b) and (c) show the 
difference between the measured and expected amplitude 
and phase envelopes for noiseless simulated data, whilst 
panels (a) and (f) show the equivalent quantities for the 
noisy data. The dotted vertical lines mark the start and 
end of the simulated null. 
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Fig. 6. Complex longitude-dependent modulation enve- 
lope and residuals for PSR B0809-|-74-like simulated data. 
See caption for Fig. |5| for descriptions of the data shown 
in each panel. 
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Fig. 8. Two-dimensional fluctuation power spectrum of 
simulated PSR B0943+10-like data. See caption of fig. |. \^ 
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Fig. 9. Complex time-dependent modulation envelope 
and residuals for PSR B0943-|-10-like simulated data. See 
caption for Fig. ^. 
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Fig. 10. Complex longitude-dependent modulation enve- 
lope and residuals for PSR B0943+10-like simulated data. 
See caption for Fig. ||. 
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Fig. 11. Power spectrum of the time-dependent ampli- 
tude envelope estimated from the noisy data. Values are 
normalised by the mean of noise-floor bins. 




-0.5 -0.49 -0.48 -0.47 -0.46 -0.45 -0.44 -0.43 
Frequency (pulses"^) 

Fig. 12. Portion of "stacked" 2DFS around the drifting 
component. Values are normalised by the mean of noise- 
floor bins. The central feature peaks at a value of ~ 487. 



